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Abstract 

It is known that classical Dirichlet-Neumann loosely coupled partitioned 
schemes for fluid-structure interaction (FSI) problems are unconditionally 
unstable for certain combinations of physical and geometric parameters that 
are relevant in hemodynamics. It was shown by Causin, Gerbeau, and No- 
bile [T^ that these instabilities are associated with the "added-mass effect" . 
The "added mass effect" was studied in ^7] by considering the simplest FSI 
problem which captures the main features responsible for the instabilities 
of classical Dirichlet-Neumann loosely coupled schemes. By considering the 
same test problem, the present work shows that a novel, partitioned, loosely 
coupled scheme, recently introduced in [TI], called the kinematically coupled 
/3-scheme, does not suffer from the added mass effect for any (3 G [0, 1], and 
is unconditionally stable for all the parameter values in the problem. Nu- 
merical results are presented for a fully nonlinear benchmark FSI problem in 
hemodynamics of Formaggia et al. [30J . The results show that the scheme is 
stable for this benchmark problem even for the parameter values well within 
the parameter range for which the classical Dirichlet-Neumann schemes are 
unconditionally unstable. The main contribution of this work is in explicitly 
showing how the implicit enforcement of the kinematic coupling condition 
and the inclusion of the structure inertia into the fluid sub-problem avoid 
the presence of the added mass effect in the kinematically-coupled /3-scheme. 
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1. Introduction 

Fluid-structure interaction (FSI) problems have important applications 
in various areas including bio-fluids and aero-elasticity. They have been 
extensively studied from the numerical, as well as analytical point of view 

[a El [3 lini [111 [111 EQl E21 Ea EZl ESI ED Ea 

A set of popular numerical schemes for FSI in blood flow includes partitioned 
schemes (loosely or strongly coupled). Partitioned schemes typically solve an 
underlying multi-physics problem by splitting the problem into sub-problems 
determined by the different physics in the coupled problem. In particular, in 
fluid-structure interaction problems the fluid dynamics and structure elasto- 
dynamics are often solved using separate solvers. In loosely coupled schemes 
only one iteration between the fluid and structure sub-problem is performed 
in each time step, while in strongly coupled schemes several sub-iterations 
between the fluid and structure sub-problems need to be performed at each 
time step to achieve stability. 

The main advantages of loosely coupled partitioned schemes are mod- 
ularity, simple implementation, and low computational costs. However, in 
P,7J it was proved that for certain combinations of physical and geometric 
parameters (which are realistic in blood flow) "classical" Dirichlet-Neumann 
loosely coupled schemes are unconditionally unstable. In the same paper, 
the authors showed that this instability is due to the "added-mass effect". 
Namely, it was shown that a portion of the fluid load to the structure in 
the coupled FSI problems can be written as an additional inertia term in 
the structure equation coming from the fluid mass (added mass). In numer- 
ical schemes in which this term appears explicitly, which is the case in the 
classical Dirichlet-Neumann loosely coupled partitioned schemes, the added 
mass term acts as a source of instabilities when the structure is too "light" 
to counter-balance the kinetic energy of the "heavy" fluid load. 

To get around these difficulties, several different loosely coupled algo- 
rithms have been proposed. The method proposed in [5] uses a simple mem- 
brane model for the structure which can be easily embedded into the fluid 
problem where it appears as a generalized Robin boundary condition. In this 
way the original problem reduces to a sequence of fluid problems with a gen- 
eralized Robin boundary condition which can be solved using only the fluid 
solver. A similar approach was proposed in [13] where the fluid and structure 
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were split in the classical way, but the fluid and structure sub-problems were 
linked via novel transmission (coupling) conditions that improve the con- 
vergence rate. Namely, a linear combination of the dynamic and kinematic 
interface conditions was used to artificially redistribute the fluid stress on 
the interface, thereby avoiding the difficulty associated with the added mass 
effect. 

A different approach to stabilization of loosely coupled (explicit) schemes 
was proposed in |15]. There, a stabilization based on Nitsche's method [M] 
was used with a time penalty term giving an L^-control on the fluid force vari- 
ations at the interface. We further mention the scheme proposed in [0] where 
a Robin-Robin type preconditioner was combined with Krylov iterations for 
a solution of an interface system. 

For completeness, we also mention several semi-implicit FSI schemes. The 
schemes proposed in |26l [H |2] separate the computation of fluid velocity 
from the coupled pressure-structure velocity system, thereby reducing the 
computational costs. Similar schemes, derived from algebraic splitting, were 
proposed in [3l HTj. We also mention [H] where an optimization problem is 
solved at each time-step to achieve continuity of stresses and continuity of 
velocity at the interface. 

Recently, a novel loosely coupled partitioned scheme, called the "kinemat- 
ically coupled /3-scheme", was introduced in [TT]. This scheme successfully 
deals with problems associated with the added mass effect in a way differ- 
ent from those reported above. The kinematically coupled /3-scheme is a 
modification of the kinematically coupled scheme first introduced in |33] . 
The parameter /3 was introduced in [11] to increase the accuracy. This pa- 
rameter distributes the fluid pressure between the fluid and the structure 
sub-problems. For /3 = the entire pressure is considered in the fluid sub- 
problem giving no pressure loading in the structure sub-problem, while for 
(3 = 1, the entire fluid pressure is used to load the structure sub-problem. 
The case (3 = corresponds to the classical kinematically coupled scheme 

The first version of the scheme, presented in [TT], concerns FSI problems in 
which the structure is described by a thin structure model (a thin membrane 
or shell model; elastic or viscoelastic). Extensions of the scheme have recently 
been developed by the authors to handle FSI problems with thick structures 
modeled by the equations of 2D or 3D elasticity [13], and to handle FSI 
problems with multiple layered structures (thin-thick-thin-thick. . . ; elastic or 
viscoelastic) [14] . Modeling FSI with multi-layered structures is particularly 
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useful in the blood flow application since the walls of major arteries are 
composed of several layers: the tunica intima, the tunica media, and the 
tunica adventitia, separated by the thin elastic laminae. 

Recently, the authors have used the kinematically-coupled scheme to 
prove the existence of a weak solution to a fully nonlinear FSI problem be- 
tween an incompressible, viscous fluid and a thin structure modeled by either 
the elastic or viscoelastic shell equations |10] . The existence proof is based on 
constructing approximate solutions using the kinematically-coupled scheme, 
and showing that the approximate solutions converge to a weak solution as 
the time-discretization tends to zero. This existence result relies on energy 
estimates which show that in the kinematically-coupled scheme the energy 
of the coupled FSI problem is well-approximated by the discretized prob- 
lem. A consequence of this result is that the kinematically-coupled scheme 
is stable. The stability of the kinematically coupled scheme is achieved by a 
clever splitting of the structure problem into the inertia part (combined with 
the viscoelastic part when the (thin) structure is viscoelastic) and the elas- 
tic part. The structure inertia is then used as a boundary condition in the 
fluid sub-problem, while the elastodynamics part is solved separately. The 
inclusion of the structure inertia (and possibly the viscoelastic part of a thin 
structure) into the fluid sub-probem was enabled by an implicit enforcement 
of the kinematic coupling condition, which enforces the no-slip condition at 
the fluid-structure interface. 

The present manuscript reveals, explicitly, how and why the kinematically- 
coupled /3-scheme does not suffer from the added mass effect for any < 
P < 1. The method of proof is different from the one presented in |50]. It 
is based on the approach presented in [T7] where it was shown that the clas- 
sical Dirichlet-Neumann loosely-coupled scheme is unconditionally unstable. 
The result in [TTj was proved on the simplest FSI problem which still re- 
tains all the main features responsible for the stability issues associated with 
the classical loosely coupled schemes. The present manuscript considers the 
same test problem for which it is shown that the kinematically-coupled /3- 
scheme is unconditionally stable for all < /3 < 1. The theoretical results are 
conflrmed by numerical simulations which where used to solve a fully non- 
linear benchmark problem in FSI in hemodynamics, originally presented in 
[30] and used for testing of various numerical schemes O HJl |3l HHl |33]. The 
hemodynamics and structure parameters in this benchmark problem fall into 
the critical regime for which the classical loosely-coupled Dirichlet-Neumann 
schemes are unstable. Our results show that the kinematically-coupled /3- 
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scheme is stable for those parameters, as predicted by the theory. We also 
show that the kinematically-couped /3-scheme compares well with the simula- 
tions obtained using a monolithic scheme by Badia, Quaini, and Quarteroni 
[3l |l6] . Further numerical results obtained using the kinematically-coupled 
/3-scheme and a comparison with experimental measurements can be found 
in [ig. 

We conclude this section by mentioning the most recent generalization of 
the kinematically-coupled scheme by Fernandez et al., called "the incremen- 
tal displacement-correction scheme" [211 12H1 [25]. In [25] the authors prove 
convergence of this numerical scheme by analyzing a simplified linear model 
problem. The incremental displacement-correction scheme threats the struc- 
ture displacement explicitly in the fluid sub-step and then corrects it in the 
structure sub-step. This scheme can also be viewed as a kinematic pertur- 
bation of a semi-implicit scheme. This observation was crucial for proving 
convergence of the scheme [25]. A different approach, however, needs to 
be used to prove stability of the kinematically coupled /3-scheme due to its 
particular splitting of the normal stress. This gives rise to the various diffi- 
culties in estimating the pressure terms which cannot be handled by using 
the classical methods presented in [25] . 

Therefore, the present manuscript provides an original look at the study 
of stability of the kinematically coupled /3-scheme. The method of proof is 
different from those used in [2S1 SD] • The results of this manuscript provide 
the information about the stability of the kinematically coupled /3-scheme 
and its relation to the added mass effect which is not captured by the results 
in [25l Bo] . This is achieved by an "explicit calculation" of the solution to the 
simplified FSI problem considered here. The main contribution of this work is 
in explicitly showing how the implicit enforcement of the kinematic coupling 
condition and the inclusion of the structure inertia into the fiuid sub-problem 
avoid the presence of the added mass effect in the kinematically-coupled {3- 
scheme. 

The manuscript is organized as follows. We begin by a description of the 
kinematically-coupled /3-scheme for a full FSI problem in which the structure 
is modeled by a cylindrical, linearly viscoelastic membrane shell model, also 
known as the ID generalized string model. See Section |2] Then, a simplified 
FSI problem is introduced, and the corresponding application of the kine- 
matically coupled /3-scheme is presented, see Section |3j An analysis showing 
the reasons why the kinematically coupled /3-scheme does not suffer from 
the added-mass effect and is unconditionally stable for any < /3 < 1, is 
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presented in Section |4j Numerical results are presented in Section [5] where a 
fully nonlinear FSI benchmark problem in hemodynamics is considered. Sta- 
bility of the scheme is shown for the parameter values well within the range 
of instability of Dirichlet-Neumann schemes. A remark discussing a com- 
parison with the classical Dirichlet-Neumann loosely coupled scheme, and a 
discussion related to the "added mass effect" are presented in Section |6| 

2. Problem definition for a fully nonlinear FSI problem 

We consider the flow of an incompressible, viscous fluid in a channel of 
radius R and length L, see Figure [TJ To fix ideas we will be assuming that 
the channel is a subset of M^, although none of the ideas related to the 
definition of the kinematically coupled /3-scheme depend on the dimension 
of the problem. The channel reference domain is denoted by f2 = (0,L) x 
{—R,R) and the lateral boundary by 

r = {{z,r) eR"^ \ < z < L, r = ±R}. 

The lateral boundary of the channel is assumed to be deformable, with a 




Figure 1: Deformed domain ^{t). 

negligible longitudinal displacement. Without loss of generality, we consider 
only the upper half of the fluid domain with a symmetry boundary condition 
at the bottom boundary r = 0. The fluid domain which depends on time is 
not known a priori. Denote by 

n{t) = {{z,r) eR'^ : z e {0,L), re (0, R + r]{z, t)} 

the fluid domain at time t, and the lateral (top) boundary by 

r(t) = {{z,r) eR^ -.r = R + r]{z, t), z E (0, L)}, 

where rj denotes the vertical (radial) displacement of the lateral boundary. 
Furthermore, let F;, = (0, L) x {0} denote the bottom portion of the boundary. 
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while Tin = {0} x (0, R), Tout = {L} x (0, R) denote the inlet and outlet parts, 
respectively. 

We consider the flow of an incompressible, viscous fluid driven by the 
inlet and outlet pressure data Pin/outit), with the fluid velocity on T{t) given 
by wBr, and the symmetry boundary condition at the bottom part of the 
boundary. Thus, the fluid problem reads: Find the fluid velocity u = 
{uz{z, r, t),Ur{z, r, t)) and pressure p = p{z, r, t) such that 

P/(^ + (u-V)u) = V-a infi(t), t e (0,T), 

V-u = infi(t), 

Ur = 0,^ = onTfeX (0,T), 

crn = -pin/out{t) on Tin/out X {0,T), 
u = wer on r(t), tG(0,T), 

where p/ is the fluid density, and cr is the fluid stress tensor. For a Newtonian 
fluid the stress tensor is given by cr = —pi + 2/iD(u), where /i is the fluid 
viscosity and D(u) = (Vu + (Vu)'^)/2 is the rate-of-strain tensor. 

The structure problem is defined by solving an elastodynamics problem 
for a cylindrical linearly viscoelastic membrane shell, also known as a ID 
generalized string model [48j, capturing only vertical (radial) displacement 
7] = r]{z,t) from the reference configuration: 

Psh^ + Co, - C,— + - ^1^^ = /' (2) 

with the boundary conditions ri{0) = r]{L) = 0, and the initial condition 
given by zero initial displacement and zero initial structure velocity. The co- 
efficients ps and h are the structure density and thickness, respectively, while 
the constants Ci,Di > 0,2 = 1,2, are the elastic and viscoelastic structural 



coefficients [IS 

The coupling between the fluid and the structure is defined by the kine- 
matic and dynamics lateral boundary conditions, respectively: 

(^(2;,t),0) = u{z,R + 7]{z,),t), (.3^1 
f{z,t) = -J{z,t){(Tn){z,R + r]{z,t),t) ■ Cr, 



where J{z, t) = —\l 1 + ^^(z, t)^ denotes the Jacobian of the transforma- 
tion from the Eulerian framework used to describe the fluid equations, to the 
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Lagrangian coordinates used in the structure equations. Here / in (|3]) is de- 
fined by the left hand-side of ([2]). The kinematic lateral boundary condition 
describes continuity of velocities at the fluid-structure interface (the no-slip 
condition), while the dynamic lateral boundary condition describes the bal- 
ance of contact forces at the fluid-structure interface: the rate of change of 
the structure momentum is a result of the balance of all the forces acting at 

2.1. The kinematically coupled ^-scheme 

To solve this problem numerically, we consider the kinematically coupled 
/5-scheme, introduced in [TT]. This scheme is based on the time-discretization 
via operator splitting known as the Lie splitting (see e.g., [32l [H])- To per- 
form the Lie splitting, problem ([l])-(|3]) must first be written in the form: 
dU/dt = AU,t e (0,T). To do this, we utilize the kinematic lateral bound- 
ary condition and express the time-derivative of rj as the trace of the fluid 
velocity on F, u|r ■ e.^. The resulting system in first-order form is given by 
the following: 

Problem 2.1. Find {u,p,ri) such that 



p^(— + (u-V)u) = V-(T mn{t), te{0,T), 



V-u = infi(t), 
Ur = 0, ^ = onr6x(0,T), 



dr 



crn 



-Pin/out{t) on Tin/out X (0,r), 

u = ^e, onr(t), t G (0,T), 



+Dou\r ■ Bj. — Di — ^ ' ■ = —Jan Br onrx(0,T). 

To deal with the motion of the fluid domain an Arbitrary Lagrangian- Eulerian 
(ALE) approach is used [Sni EH 112] • We introduce a family of (arbitrary, in- 
vertible, smooth) mappings At defined on the reference domain Q such that, 
for each t G [0,T], At maps the reference domain Q onto the current domain 
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Q{t). Now, ALE time-derivative of function / defined on Q{t) x (0,T) is 
given by 

dl 

dt 

where w denotes the domain velocity given by 



l+w.V/, (4) 



dAt 

w(2;,r, t) = —^{z,f), where (-2, r) = A^^{z,r). (5) 

Problem Q-Q can now be written in the following first-order ALE form 
(for more details please see [H]): 

Problem 2.2. Find {u,p,ri) such that 



((u-w)-V)u) = V-cr inl](t), t e (0,r), 

V-u = mQ{t), 
du 

= 0, ^ = onrfeX(0,T), 



crn 



dr 

= -Pin/out{t) on Tin/out X {0,T), 

u = ^e, onr(t), t G (0,T), 

+DQu\r ■ e.r - Dl — ^ f ■ e,. = -Jan - Br onrx(0,T). 

The strategy of the kinematically coupled scheme is to split this problem 
into a fluid sub-problem and a structure sub-problem in such a way that 
the inertia of the thin structure is coupled with the fluid sub-problem via a 
"Robin-type" boundary condition on T{t). When the structure is viscoelastic, 
the structural viscosity can be treated together with the structure inertia as 
a part of the same boundary condition (i.e., all the terms in the structure 
equation involving the trace of fluid velocity on F can be used as a boundary 
condition for the fluid sub-problem). The elastodynamics of the structure 
problem is solved separately in the structure sub-problem. 

Additionally, in the /3-scheme the fluid stress is split into two parts: 

crn = crn + j3pn —jSpn, j3 G [0, 1], 
(I) ill) 
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where the case /3 = corresponds to the classical kinematically coupled 
scheme introduced in |33] . Part I of the fluid stress is used in the fluid sub- 
problem, while Part II of the fluid stress (the /3-fraction of the pressure) is 
used to load the structure in the elastodynamics sub-problem. Thus, the 
dynamic coupling condition 

d(\i\ ) d'^ri [ 9^ful \ \ 

Psh^^-Br = -CqT] + Ci-^+ iDou\r - Di + Jan + J^pn \ ■ e, 

—J(3pn ■ Br 
is split into the fluid part 

d(u\ ) (9^ful ) 

PART 1 : psh — = — -Do^lr " + -Di — — ■ + Jan ■ + jPpn ■ e^, 

ot oz^ 

and the structure part 

PART 2 : Psh^^^ ■ e, = -CoV + C,^-J/3pn, with p,/i^^-e, = p,h^. 
The scheme is then defined by the following. 

The fluid sub-problem: 

Step 1. The Stokes problem is solved on a fixed fiuid domain determined 
from the previous time step. The boundary condition at lateral boundary 
is given by PART 1 of the dynamic boundary condition, with Ppn given ex- 
plicitly from the previous time step. The displacement of the structure stays 
intact. The problem reads: 

Given p", r/" from the previous time step, find u, p, rj such that for t G 

P;^|.=v-o-, v-u = o i22^](n 



Ps/i-^lr ■ e,. + -Dou|r ■ - -Di^|r ■ e^ + W (1 + I — I [an) ■ er 



+ P"n-er 0:2 (0,L), 

u ■ ez = on (0, L), 
^ 1=0 iMO.L). 
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with the following boundary conditions: 

du 



dr 



z, 0, t) = Ur{z, 0, t) = on Tf,, = on T{f 
u{0, R,t) = u{L,R,t) = 0, 



crn 



Pin {t)n\in on Tin, crn 



out 



-Pout{t)n\out on Tout, 



and initial conditions 



u(r) = u", r/(^) = r/^ 

Then set u^+^Z^ ^ ^(r+i), r^^+^Z^ = r7(r+i), = p(r+i). 
The advection problem: 

Step 2. Solve the fluid and ALE advection sub-problem deflned on a flxed 
domain fi(t"), with the ALE velocity w""*"^/^ deflned by ^ based on the 
domain Q{t^) and the corresponding ALE mapping. The problem reads: 
Find u and r] such that for t E (t", t""*"^) 



dfj 

di 



+ (u"+V3 _ ^"+1/3) . Vu = 0, in n{f 
on{0,L), 



du , 



Ps/is-^|r = 0, oii(0,L), 



with the following boundary conditions: 

u = u"+i/=^ Oil r^'/l where 



{xGM^|xGa^](^),(u 



n<0}, 



and the initial conditions 
Then set u"+2/3 = u(r+^), = r/(r+i). 



The structure (elastodynamics) sub-problem: 
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Step 3. Step 3 involves solving the elastodynamics problem for the location 
of the deformable boundary by involving the elastic part of the structure 
which is loaded by Part II of the normal fluid stress. Additionally, the fluid 
and structure communicate via the kinematic lateral boundary condition 
which gives the velocity of the structure in terms of the trace of the fluid 
velocity, taken initially to be the value from the previous step. The problem 
reads: 

Find u and rj, with computed in Step 1 and 77" obtained at the previous 
time step, such that for t G (t", t""*"^) 

dr] 

— = u\r(t)-er on{0,L), 



Psh^\r ■^r + CoV-C,^,= (1 + (^) V+^n ■ on (0, L), 
with the boundary conditions: 

V\z=0,L = 0, 

and the initial conditions: 

U(r) = u"+2/3^ ^(^fi^ = ^n+2/3_ 

Then set u"+i = u(r+^), r/'"+^ = r]{t''+^). 
Do = t"^^ and return to Step 1. 

A block-diagram showing the main steps of the scheme is given in Fig- 



ure 2.1 More details about the scheme can be found in ITT]. 



3. The simplified test problem 

We show that the kinematically coupled /5-scheme is unconditionally sta- 
ble when applied to a simplified problem, introduced in [T7j as the simplest 
problem which captures the main features related to the instabilities in clas- 
sical loosely-coupled schemes caused by the "added mass effect". The sim- 
plified problem consists of solving the time-dependent Stokes equations for 
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^^^^^^^^^^^^^^^^^^^^^^^^^ 

Initial structure position and 
, fluid pressure , 



IVIESH UPDATE ] 

o 



FLUID SOLVER 

witli boundary condition on T 

— + o"n"-' = -pp''-'n"-' 
dt 



STRUCTURE SOLVER 

with pressure forcing / = j8 

and initial structure velocity = trace of fluid velocity 



4 a=n+i 



Figure 2: A block diagram showing the main steps of the kinematicahy coupled /3-scheme. 



an incompressible, inviscid fluid in a 2D channel with deformable walls, and 
with the elastodynamics equations given by the ID generalized string model. 
Moreover, it is assumed that the displacement of the deformable wall is small 
enough so that it can be neglected in the fluid flow problem. In this case 
the geometry of the fluid domain is fixed while the small deformation of the 
boundary is calculated using the elastodynamics equations which are coupled 
with the fluid flow through the kinematic and dynamic coupling conditions. 
Thus, the problem is defined on the reference fluid domain Q with the lateral 
boundary F. 

We will be assuming that the inlet and outlet pressure data, which are 
functions of time, are uniformly bounded. This is a reasonable assumption for 
the blood flow application, since the pressure is typically a periodic function 
of time, with a bounded amplitude. 

The fluid problem reads: Find the fluid velocity u = u{z,r,t) and 
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pressure p 



p{z, r, t) such that 
du 

V u = 
u ■ n = 

P = 

U ■ Br = 








Pin/ outify 
W 



inQx (0,T), 

in X (0,T), 
on X (0,T), 
on Tin/out X (0,T), 
on r X (0,T), 



(6) 



with the initial velocity and pressure equal to zero. 

The structure problem is defined by solving for t] 
generalized string model: 



/, 



ri{z,t) the ID 



(7) 



with boundary conditions ?7(0) = rj{L) = 0, and initial conditions given by 
the zero initial displacement and zero initial structure velocity. 

Remark 3.1. Notice that the simplified problem does not have the viscoelastic 
terms. The addition of the viscoelastic terms only contributes to the stability 
of the scheme. For more details about the analysis of the fully nonlinear FSI 
problem with and without the viscoelastic terms see IJI^ . 

The coupling between the fiuid and structure is defined by the kinematic 
and dynamic lateral boundary conditions, which, in this simplified problem, 
read: 

/ = P|r, 



(8) 



The coupled fluid-structure interaction problem can be written as follows: 
Problem 3.1. 



Pf 



du 

dt 



+ Vp = 





V u = 





u ■ n = 





P = 


Pin/ 




drj 


U ■ Br = 


di 


P = 


Psh 



+ Cor] - C, 



mQx (0,T), 

inVlx (0,T), 
onTh X (0,T), 
on Tin/out X (0,T), 

onT X (0,r), 
onT X (0,T). 



(9) 
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P = P^h—\^-er + CoV-Ci^ onrx(0,r). 



To perform the Lie splitting, this system is written as a first-order system 
by using the kinematic couphng condition to obtain: 

Pf-Ql + ^P = infi X (0,T), 

V-u = infi X (0,r), 

u • n = on Tfe X (0,T), 

P = Pin/outif) on ri„/o„j X (0,T), 

u-Br = ^ onr X (o,r), 

du, ^ ^ d'^V 

(10) 

As before, this point is crucial in order to perform the splitting which gives 
rise to a stable scheme. Namely, in the fiuid sub-problem, which we write 
next, the structure inertia will be taken into account implicitly in the bound- 
ary condition on F, which is a crucial ingredient for the stability of the 
/3-scheme. 

Before we write the main splitting steps we also notice that the fiuid 
stress, which in this simplified example is only the pressure, will be split into 
two parts: Part I which is given by p — and Part 2 which is given by [5p 
so that p = {p — (3p) + /3p=Part 1 + Part 2. 

Thus, the main steps of the Lie splitting for the simplified FSI problem 



(10) are: 



The fluid sub-problem: 

Step 1. The Stokes problem is solved on a fixed fiuid domain, with a bound- 
ary condition which couples the structure inertia with the Part I of the fiuid 
stress. Here, the portion f3p of the stress is taken explicitly, while the rest is 
taken implicitly. The displacement of the structure stays intact. The prob- 
lem reads: Given from the previous time step, find u, 7] and p such that 

forte (r,r+i); 

du ^ ^ 

Pf-Q^ + ^P = mn, 

V ■ u = inVt, 
u ■ n = onTb, (11) 

P = Pin/out(t) on Tin/out, 

du I 

p-fip^ = psh—\^-er onT. 
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The structure displacement stays intact, namely 

1=0 on r. (12) 
The initial conditions are taken from the previous time step. Then set 

u"+V2 ^ = u|r(r+^),p"+' (13) 



The structure sub-problem: 

Step 2. The structure problem is solved with the pressure load Pp^^^ 
just calculated from the fluid sub-problem: Find u and f] such that for 
t e the following holds: 

Psh—\^-er + Cor)-Ci^^Pp^+'\r on{0,L), 

— {z,t)=u\r-e., on{0,L) 

The Euid velocity in fl remains intact, namely: 

du ^ ^ 
— = in Q. 

dt 

The initial data correspond to the solution of the problem in Step 1 evaluated 
at Then set 



The pressure formulation: 

In this simplified model, the problem is Step 1 can be entirely formulated 
in terms of the pressure and the normal trace of the fluid velocity on F. 
Namely, by taking the divergence free condition in the differentiated first 
(momentum) equation, the problem can be written in terms of the Laplace 
operator for the pressure. Similarly, we can re-write the boundary conditions 
on r in terms of the pressure and the normal trace of the velocity on F. 

To simplify notation we will be using n to denote the normal, e^, on F, 
and Ur\r to denote the normal trace of the fluid velocity on F: 

n :— er on F and Urlv '■— u|r • e^. 
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The fluid sub-problem in terms of p and Ur\r 
Step 1'. Find p and Ur\r such that for t G (t", t"+-^); 



p + 



-Ap 
p 
dp 

psh dp 
Pf dn 



0, 

Pin/out(j^^ 





in Q, 

O-fi '^in/outi 

on Tb, 

on r, 



du. 



T r 



1 dp\i 



dt 



on r, 



Pf dn 

The first four equations in Step 1', supplemented with the initial data 
for the pressure determine a time-dependent Robin problem for the pres- 
sure. The time-dependence enters through the inlet and outlet pressure data, 
which are functions of time. Notice that this problem already incorporates 
a portion of the coupling of the underlying FSl problem. This problem 
has a unique solution in the space C{0,oo; H^{Q)) provided that the data 
Pin/out e C(0, oo; Hy\T,^/out)), and p" G C(0, oo; R-y^T)). 

One can also show (see e.g., [T7j) that the problem for the structure, 
defined in Step 2, has a unique solution r] G C((0, oo); ffg (0, L)) provided 
that the initial data are such that ri\t=o G -f^o(0, L), and drj/dt\t=Q G -f^o(0, L), 
with p"+i|r G C(0, oo; H^/^{T)). 

To simplify notation, as in [17], we introduce a linear, symmetric, posi- 
tive definite operator £ defined by the elasticity tensor associated with the 
structure problem. Namely, for z/,^ G Hq{0,L), we define 



< Cr],^ >:= as{ri,^), 
where as{ri,C,) is the inner product on H^{0,L) defined by 



(15) 



CoTj^dx 



V ^^^^^ 







' dx dx 

With this notation, the structure equation in Step 2 can be written as 

d'^rj 



Psh 



dt^ 



Ct] = I3p\r. 



Operator C can be extended to the full structure operator defining the Koiter 
shell model, described in ITGl. 
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We now summarize the main steps of the sphtting sheme. With a shght 
abuse of notation, we emphasize the discretized form of the Lie sphtting, 
defining our kinematicaUy coupled /3-scheme: 

Step 1. Given <|r, p",??"*, Gnd p"+V2^ Ur^^^^\r, 77^+^2, such that for t e 



p 



pn+l/2 
Qpn+l/2 

dn , 

n+1/2 I Psh dp^-^y^ 
PS dn 



— Pin/outit) 

= 



in Q, 

on Tft, 
onV, 



dt 



,rt+l/2 



1 ap"+V^|r 



012 r, 
on r. 



^ i / = Tj 

with the initial data given by u^\r, p", rj^- 

Step 2. Given p''+^/^ , Ur^^^^\r, r]''+^'^ find p''+\ <+V, such that for 

t e (r,r+i): 

a2 n+l 



drj 



.n+l 



dt 



on r 

on 



witii the initial data given byp^'^^^'^, Ur'^^^'^\r, Ty^+^Z^ ^n+i_ 

We study the stability of the kinematicaUy coupled ^-scheme for this 
problem next. 



4. Stability analysis 

To study the stabihty we introduce an operator V : iJ"^/^(r) — > Q, where 

Q = {qe H\Q) I g|r.„/„„^ = 0}, (16) 
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such that 



-AVw 
Vw 
dVw 



Vw + 



dn 

Psh dVw 



0, in Q, 

on Tin/ out ^ 

on Ffe, 

w on r. 



(17) 



Pf dn 

Thus, operator V associates to every w E H~^/'^(T) the solution of the 
pressure problem in Step 1', with the homogeneous inlet and outlet data 

Pin/out 0. 

We are interested in the trace on F of this pressure solution. For this 
purpose, introduce the operator S : H^^^'^{T) — )• H^^'^{T) by 



(18) 



Sw = Vw\ 

One can prove the following standard results (see [IZl): 
Proposition 4.1. The operators satisfies the following properties: 

1. The operator S : H^^^'^{T) — t- H^^'^{T) is continuous. 

2. The operator S is compact, self-adjoint, and positive on L'^{T). 

Proof. Since this result is quite standard, we just outline the proof. We start 
by writing the weak formulation of problem (17) which reads: find Vw G Q 
such that 



VVw : V0 + 



PL 
Psh 



Vwc 



Pl_ 
Psh 



{w,(t)), V0 G Q, 



(19) 



where {w, (p) is a duality pairing between H-^/^{T) and H^/^(T). The first 
assertion (continuity) follows directly from here. 

To show the second statement of the proposition, we consider w,v E L'^{T) 
and notice: 



Swv 



Vwv 



psh 
Pf Jn 



Wv : Ww 



VvVw. 



(20) 



Here, the first equality follows from the definition of operator S. For the 
second equality we notice that Vw G Q is an admissible test function and, 
therefore, the second equality is just formula (19) with the test function (j) = 
Vw. Now, we consider V as an operator on L%F), i.e., S : ^^(F) ^ ^^(F). 
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For every w G L'^iV)^ we have Vw G Q. Now, the trace theorem imphes 
that Sw = Vw\r G H^/^{T) and so we have lm{S) C H^/^{T). Since T 
is a bounded set, embedding if^/^(r) ^ -^^(r) is compact, and thus S is 
a compact operator. Formula (20) imphes that operator S is positive and 
self-adjoint, which completes the proof. □ 
To study the solution of the corresponding non-homogeneous problem 
with p = Pin/out(t) on Tin/ out, let us introduce an arbitrary continuous ex- 



tension operator Ep : H^/'^{dn\T) H^{n) so that Epq] 



dn\r 



q and 



\EFq\\m(n) < C'||'?||_H'i/2(ar2\r)- Let p* G C{0, oo; H^{Q)) be the solution to 



p + 



-Ap* 


= AEfP, 


p 


= 


dp* 


dEpp 


, Qn 




Psh dp* 


d'^p 


Pf dn 


dn 



in Q, 

on r 



in/ out 1 



on Ffe, 
on F, 



where p = Pin/out on Tin/out- Now, the solution to the pressure problem in 
Step 1' is given by 

p = p* + EFp + Vi(3p''). 

Let us denote 

Pext =P*\r + Epplr- 
Now we can write the trace of the pressure solution in Step 1' on F as 

p|r=Pext + 5(/3p"). (21) 

This trace of the pressure is used to load the equation for the structure in 
Step 2. Thus, the structure problem in Step 2 can now be written as: Find 
7] such that 

Psh^ + Cr^ = mtt + S{m)- (22) 



One way to discretize this equation in time is to use an implicit marching 
scheme 



pji 



n+l 



2^" _|_ 



At2 



+ Cr] 



n+l 



/3(p:+i+5(/3p")). 



(23) 



This scheme is unconditionally stable (provided that the right hand-side con- 
verges) . 
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We chose to discretize equation (22) in time using a 6'-scheme discussed 
in [32]: 



Psh 



n+l 



2r/" + r] 



n-l 



At2 



+ £(077"+i + (1 - 2^)77" + Or] 



"'') = /3(pLY + 5(/3p")). 

(24) 

It was shown in [32] that, for a given fixed right hand-side (source term), 
this scheme is stable for all < ^ < 1/2. Thus, we have unconditional sta- 
bility with respect to the arbitrary ratios of the fluid and structure densities, 
provided that the right hand-side of this equation converges as n — > 00. 

A crucial point to observe here is that the right hand-side of this equation, 
which comes from the pressure loading, is given by an iterative procedure, and 
can be written entirely in terms of the initial pressure, the external pressure 
(Pin/out), and the operator S whose maximum eigenvalue, as we shall show 
below, is always less than 1, for all the choices of pf and pgh. Moreover, 
we will show below that the right hand-side converges, as the number of 
iterations n— )'OoifO</3<l, for all the choices of pf and psh. 

To analyze the right hand-side of equation (22), we first study the eigen- 



values of the operator S. As we shall see below, it is convenient to express 
the eigenvalues of S via the eigenvalues of the "Neumann to Dirichlet" oper- 
ator Ma '■ H~^^'^(T) — !■ i7^/2(r) which is defined to be the trace on T of the 
operator 7^ : R-^/^T) Q: 



where 71 associates to every w 
following (pressure) problem: 



r -AUw 
TZw 
dTZw 



Maw = nw\r, (25) 
G H~^/^{T) the solution TZw e Q of the 



dn 
dTZw 

dn 



0, in Q, 

on Tin/outi 

on Tb, 

w on r. 



It can be shown (see [T7]) that operator Ma '■ H^^^'^{T) — )■ H^^'^{T) is contin- 
uous, and that Ma is compact, self-adjoint, and positive on L'^{T). Moreover, 
the eigenvalues pi of Ma are decreasing to zero {pi = L/{mth.{iT[R/L)), i = 
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1,2,...), with the maximum eigenvalue fimax given by 



vrth 



/tcR 



We will use this knowledge to calculate the eigenvalues of operator S. Let 
/i be an eigenvalue of operator Ma- Then there exists a vector v such 
that 

Mav = fJ^v- 



Recall that, by definition of J^a, M.av = T^f |r and 



we calculate 



dn 



V. Using this, 



nv\r + 



Psh dTZv I 



(/i H )f r- 

Pf 



Pf dn 

This implies that TZv also satisfies the following Robin problem: 



-ATZv 
TZv 
dUv 

dn 
Psh dTZv 

Pf dn 



0, 






in Q, 
on r 



in/ out 1 



on Tb, 



\p A )v on 1 . 

Pf 



But this is precisely the problem defined by the operator V in (17), with the 
data on V given by (/i + ^)v. Thus: 

Tlv = V(p+ —)v, 
Pf 

and therefore, the traces on T satisfy: 



Mav = S(p+ =(p + ^)Sv. 

Pf Pf 



Since J^av = pv we finally get 



f , P^^\ c 
pv = [p -\ )OV. 

Pf 
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Therefore, v is also an eigenvector for S, and the corresponding eigenvalue A 
satisfies: 

A = — (26) 

^ PS 

Thus, we have shown that the eigenvalues Aj of S can be expressed using the 
eigenvalues \iioi M.j^ as 

A^i H 

Pf 



We now use this information to study the right hand-side of equation (24 ). 
Since S is compact, there exists an orthonormal basis of i^^(r) composed of 
the eigenvectors {^Zj\ of S. We thus expand the solution r] and the external 
pressure data "Pext and p° in this basis: 



Then, from (24), for each j, the Fourier coefficients satisfy the following 
equation: 

P./^^^^^^^^^^^^l^^^^^ + mv^'), + (1 - 29)ivnj + Oiv"-'),) 

(28) 

The right hand-side of this equation is equal to 

mit% + ELi /3^-^^a;.(p:,v-^- + r+^Af ^(po),. 

As 77, — )■ oo, the series that defines the right hand-side converges if 

|/3A,|<1. 



From (27) we see that all Xj are strictly less than one, which implies that the 
right hand-side converges if 

< /3 < 1. (29) 
We have shown the following result: 

Theorem 4.1. The kinematically coupled ^-scheme applied to a class of FSI 
problems represented by the simple benchmark problem (10), is uncondition- 
ally stable for each /3 such that < (3 < 1. 
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5. Numerical results 

We present two numerical examples that show stability of the kinemat- 
ically coupled /3-scheme for a fully nonlinear FSl problem. The parameter 
values in the two examples are well within the range for which the classi- 
cal Dirichlet-Neumann schemes are unstable. Numerical simulations in both 
cases are preformed on a benchmark problem by Formaggia et al. [30j used 
for testing of several FSI algorithms in hemodynamics applications [5], HSj |3l 
SHI [33]. 

It was shown in fTT] that the classical loosely coupled, Dirichlet-Neumann 
scheme is unconditionally unstable if 

Pf^^max 

where Pmax is the maximum eigenvalue of the added mass operator given by 

L 

^"^'^^ " 7rtanh(f )■ 

For the parameters given in Table [1} the value of pmax is 7.46, so the critical 
value for the structure density is ps = 74.6 g/cm^. Therefore, the classical 
Dirichlet-Neumann scheme is unconditionally unstable if ps < 74.6 g/cm^. 
In Example 1, we take the density of the structure to be 70-times smaller, 
given by the physiologically relevant value of ps = 1.1 g/cm^. We compare 
our results with the results obtained using a monolithic scheme by Badia, 
Quaini and Quarteroni [Sj US] showing excellent agreement. In Example 2, 
we choose an even smaller structural density, ps = 0.55 g/cm^, and show the 
convergence of our scheme in time, confirming the stability of the kinemati- 
cally coupled /3-scheme even for the fully nonlinear FSI problems. 

The benchmark problem [30] : 

The benchmark problem consists of solving the fully nonlinear FSI prob- 
lem ([l])-(|3]) with the values of the coefficients for the fluid problem given in 
Table [1} and the values of the structural coefficients given in Table [2j 



(30) 



The flow is driven by the time-dependent pressure data: 
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Parameters 


Values 


Parameters 


Values 


Radius R (cm) 


0.5 


Length L (cm) 


6 


Fluid density pf (g/cm'^) 


1 


Dyn. viscosity /i (poise) 


0.035 


Young's mod. £'(dynes/cm^) 


0.75 X 10^ 


Wall thickness hg (cm) 


0.1 


Poisson's ratio a 


0.5 





Table 1: Geometric parameters, and fluid and structural parameters that are used in 
Examples 1 and 2 presented in this section. 



CoefRcients Values 


CoefRcients Values 


Co 4 X 10^ 


Ci 2.5 X 10^ 


Do 


Di 0.01 



Table 2: The values of the coefficients in ^ used in Examples 1 and 2 below. 



where Pmax = 2 x lO'' (dynes/cm^) and tmax = 0.005 (s). The graph of the 
inlet pressure data versus time is shown in Figure [3j The inlet and outlet 



X 10 




0.002 



0.004 0.006 0.008 0.01 

time [s] 



0.012 



Figure 3: The inlet pressure pulse for Examples 1 and 2. The outlet pressure is kept at 0. 



boundary conditions are the absorbing boundary conditions: 



at z = 




at z 



(32) 



(33) 



To deal with the motion of the fluid domain, the ALE mapping we used in 
the simulations was defined to be the harmonic extension of the mapping 
that maps the boundary of Q to the boundary of Q{t) for a given time t. 



Parameter (3, introduced in (29), which appears in Step 1 and Step 3 of 
our numerical scheme was taken to be /3 = 1. It was numerically observed 
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in [TT] that the change in /3 is associated with the change in accuracy of the 
scheme (not the stabihty), where the value of /3 = 1 provides the highest 
accuracy for this benchmark problem. We believe that the main reason for 
the gain in accuracy at /3 = 1 is the strong coupling between the fluid pressure 
(which incorporates the leading effect of the fluid loading onto the structure) 
and the structure elastodynamics, which is established for /3 = 1 in Step 3 of 
the splitting, described above. 

5.1. Example 1: ps = 1.1 g/crn?. 

We solved the fully nonlinear benchmark problem described above over 
the time interval [0,0.012]s as in [SO]- The numerical results obtained us- 
ing the kinematically coupled /3-scheme were compared with the numerical 
results obtained using the monolithic scheme by Badia, Quaini and Quar- 
teroni [3|,|l6]. Figures |4| |5] and [6] show a comparison between tube diameter, 
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Figure 4: Example 1: Diameter of the tube computed with the kinematically coupled 
/3-scheme (solid line) and with the monolithic scheme used by Quaini in [31 ES] (dashed 
line). The time step At = 10~^ is used in both cases. 

fiowrate and mean pressure, respectively, at six different times. These results 
were obtained on the same mesh as the one used for the monolithic scheme 
in [16], containing 31 x 11 Pi fluid nodes. More preciesely, we used an isopara- 
metric version of the Bercovier-Pironneau element spaces, also known as the 
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Figure 5: Example 1: Flowrate computed witli the Itinematically coupled /3-scheme (solid 
line) and with the monolithic scheme used by Quaini in [31 |3S] (dashed line) . The time 
step At = 10^^ in both cases. 
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Figure 6: Example 1: Mean pressure computed with the kinematically coupled /3-scheme 
(solid line) and with the monolithic scheme used by Quaini in [3l|46] (dashed line). The 
time step At = 10^^ is used in both cases. 

Pi-iso-P2 approximation in which a coarse mesh is used for the pressure (mesh 
size hp) and a fine mesh for the velocity (mesh step hy = hp/2). 
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The time step used was At = 10~^ which is the same as the time step 
used for the monohthic scheme. Due to the sphtting error, it is well-known 
that classical splitting schemes usually require smaller time step to achieve 
the accuracy comparable to that of monolithic schemes. However, the new 
splitting with /3 = 1 allows us to use the same time step as in the monolithic 
method, obtaining comparable accuracy, as shown in [TT] and in Figure [Tj 

The reference solution was defined to be the one obtained with At = 10~^. 
We calculated the relative and absolute errors for the velocity, pressure 
and displacement between the reference solution and the solutions obtained 
using At = 5 X 10~^, 10~^, 5 x 10~^ and 10~^. Table [s] shows the relative er- 
ror and the convergence rates for the pressure, velocity, and displacement ob- 
tained by the kinematically-coupled /3-scheme. The graphs in Figure [7] show 
the absolute error and the convergence in time of the kinematically-coupled 
/3-scheme and the monolithic scheme by Badia, Quaini and Quarteroni [3l HB] . 



At 


|p-Pre/ 11^,2 


order 




order 


l»7-'7re/lli2 


order 


|Pre/lli,2 


ll"re/lli,2 


ll'7rc/lli2 


10-" 


0.0251 




0.0223 




0.0392 




5 X IQ-'' 


0.013 


1.35 


0.0151 


0.56 


0.0175 


1.1 


lO-'' 


0.0024 


1.04 


0.0038 


0.87 


0.0038 


0.92 


5 X 10-'' 


0.0011 


1.14 


0.0017 


1.12 


0.0017 


1.13 



Table 3: Example 1: Convergence in time calculated at i = 10 ms. 



5.2. Example 2: Ps = 0.55 g/crn? . 

We consider the same test case as in Example 1, but now with the value 
for the structure density ps = 0.55 g/cm^. The problem is solved on the time 
interval [0,0.012]s and with the same mesh as the one used in Example 1. 



Figures [8} [9] and [10] show the values of the tube diameter, flowrate and mean 
pressure, respectively, at six different times, obtained with At = 10~^. 

Finally, to show the convergence of the scheme, we define the reference 
solution to be the one obtained with At = 10~^. We calculated the relative 

errors for the velocity, pressure and displacement between the reference 
solution and the solutions obtained using At = 5 x 10~^, 10~^, 5 x 10~^ 
and 10~^. Table |4] shows the convergence rates for the pressure, velocity, 
and displacement obtained by the kinematically-coupled /3-scheme. Thus, 
we see that the scheme performs well for a range of parameters for which 
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Figure 7: Example 1: The figures show absolute errors compared with the monolithic 
scheme. Top left: Error for fluid velocity at t=10 ms. Top right: Error for fluid pres- 
sure at t=10 ms. Bottom: Error for displacement at t=10 ms, |TT]. Both schemes are 
first-order accurate in time, with comparable accuracy (notice the higher accuracy of the 
kinematically-coupled /3 scheme for the displacement). 



the classical loosely-coupled, Dirichlet-Neumann scheme is unconditionally 
unstable. 



At 


|p-Pre/lli,2 


L-^ order 


1"— "rc/lli,2 


order 


l^-'7re/lli,2 


order 


|Prc/lli2 


ll"re/lli,2 


l^re/lli,2 


10~* 


0.0239 




0.0427 




0.0749 




5 X 10"-'' 


0.0096 


1.32 


0.0286 


0.58 


0.0408 


0.87 




0.0017 


1.06 


0.0067 


0.9 


0.0079 


1.021 


5 X 10-*^ 


7.72G - 04 


1.15 


0.0031 


1.13 


0.0035 


1.16 



Table 4: Example 2: Convergence in time calculated at t = 10 ms. 
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Figure 8: Example 2: Diameter of the tube computed with the liinematically coupled 
^-scheme, obtained with the time step At = lO^"* . 
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Figure 9: Example 2: Flowrate computed with the kinematically coupled /3-scheme, ob- 
tained with the time step At = 10"''. 



6. Conclusions 

We conclude this manuscript by comparing the kinematically-coupled (3- 



scheme, summarized in (22), with the classical Dirichlet-Neumann scheme 
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Figure 10: Example 2: Mean pressure computed with the Icinematically coupled /3-scheme, 
obtained with the time step At = 10^*. 



exhibiting the "added mass effect." Equation (22 ) can be solved in an imphcit 



way (written with a shght abuse of notation) as: 

a2 n+l 

Psh^^ + l^rT^' = mtl + S{m)- (34) 

The corresponding equation resulting from the classical Dirichlet-Neumann 
scheme can be written as follows. We recall that the Dirichlet-Neumann 
scheme solves the FSI problem ([6])-([8]) by solving the fluid sub-problem ([6]) 
with Dirichlet boundary data for the fluid velocity on F, given in terms of 
the structure velocity drj/dt, calculated from the previous time step, and 
then uses the fluid stress, calculated in problem ([6]), to load the structure 
in sub-problem ([t]). Figure [6] shows a block diagram summarizing the main 
steps. 

Using similar ideas as already presented in this manuscript, it was shown 
in [17] that this kind of partitioned approach leads to solving the problem 
for 7] of the form 

P^^-Q^ + ^V = P\r, (35) 

where 

p\r = Pext - P/Ma-q^- (36) 
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Initial structure position and 
structure yelocity 



MESH UPDATE 



FLJID SOLVER 



with Diriehlet boundary condition on r 



U"\ =■ 



STRUCTURE SOLVER 

with entire normal fluid stress as forcing 



- | n = n+l I 



Figure 11: A block diagram showing the main steps of a Dirichlet-Neumann scheme. 



Here p^xt comes from the "external" pressure data (pm/out), and J^a is the 



'added mass operator" introduced in the previous section by formula (25). 



While in the Dirichlet-Neumann loosely coupled scheme the w in the 



definition of operator TZ (formula (26)) corresponds to the structure inertia 
d'^rj/dt'^, in our kinematically coupled /3-scheme w corresponds to the fraction 
of the fluid stress determined by Pp. When f3 = 0, which corresponds to 
the original kinematically coupled scheme, this term is zero. We emphasize 
again that the structure inertia in the kinematically coupled schemes is taken 
implicitly in the Robin boundary condition, and not explicitely, as in the 
classical loosely coupled schemes. 



With the added mass operator the problem for the structure (35) can 
now be written as 



Pext- 



(37) 



Since in the Dirichlet-Neumann loosely coupled partitioned schemes the pres- 
sure is calculated in the fluid sub-problem using the structure velocity from 



the previous time step, this implies that, with a slight abuse of notation, (37) 
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can be written as 



q2 n+l d'^Tl^ 



More precisely, equation (38) means 



„n+l _ n n _i_ „n-l „n _ n n-1 _i_ „n-2 

7 7 "^'1 ' 'I I . ^ 7 ^7 "T 7 I n /on\ 

(Atp +P/-^A ^^^p + Cr] =pext- (39) 

It was shown in [1^ that this scheme is unconditionally unstable if 

Psh 



PfPr. 



< 1, (40) 



where pmax is the maximum eigenvalue of the added mass operator J^a- It 
was also shown in [T7j that the maximum eigenvalue Pmax is associated with 
the aspect ratio of the fluid domain R/L. The smaller that aspect ratio (the 
more slender the domain), the larger the maximum eigenvalue pmax- 

Going back to equation ( p4| ), one can now see that there is no added mass 
operator associated with the kinematically coupled /3-scheme. The right hand 



side of equation (34) is an iterative procedure which can be written in terms 



of the initial pressure and of the external (in/out) pressure data. There is 



no added mass operator Ma that appears in (34). Thus, for < /3 < 1, the 
kinematically coupled /3-scheme does not suffer from the added mass effect 
for any of the choices of psh, pf, and domain sizes. 
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